#ifndef _BML_FAST_MATH_H_
#define _BML_FAST_MATH_H_

/*! \addtogroup fast_math Fast 32-bit maths
 * @{
 * \brief Fast conversions.
 *
 * The VM creates its own environment, which is used to (un)serialize programs.
 * The map supports indexed access to keys and values and random access by key lookup.
 */

/* 32-bit floats and ints only ! */
/*! \brief Magic converting bias, hex version. */
#define BIAS_HEX ( ((23+127)<<23) + (1<<22) )
/*! \brief Converting bias as an int */
#define BIAS_INT ((long int)BIAS_HEX)
/*! \brief Converting bias as a float */
#define BIAS_FLOAT ((float)12582912.0f)

/*! \brief union to reinterpret bits in a word as an integer and a float. */
union _intfloat_conv {
	long int i;	/*!< \brief bits as integer */
	float f;	/*!< \brief bits as floating point */
};

/*! _brief Public type for conversion union. */
typedef union _intfloat_conv _IFC;

/*! \brief Convert from float to int. */
static inline long int f2i(float f) {
	_IFC c;
	c.f=f+BIAS_FLOAT;
	return c.i-BIAS_INT;
}

/*! \brief Convert from int to float. */
static inline float i2f(long int i) {
	_IFC c;
	c.i=i+BIAS_INT;
	return c.f-BIAS_FLOAT;
}


/*@}*/

#include <math.h>

#ifndef M_PI
#define M_PI ((float)3.1415926535897932384626433832795028841971693993751f)
#endif

/* adapted from http://www.devmaster.net/forums/showthread.php?t=5784 */
static inline float sine(float x)
{
    #define B ((float)(4.f/M_PI))
    #define C ((float)(-4.f/(M_PI*M_PI)))
    /*const float B = 4/pi;*/
    /*const float C = -4/(pi*pi);*/

    float y = B * x + C * x * fabs(x);

    #ifdef EXTRA_PRECISION
    //  const float Q = 0.775;
        const float P = 0.225;

        y = P * (y * fabs(y) - y) + y;   // Q * y + P * y * abs(y)
    #endif
    #undef B
    #undef C
    return y;
}



static inline float cosine(float x) {
	x+=M_PI/2;
	if(x>M_PI) {
		x-=2*M_PI;
	}
	return sine(x);
}

/* adapted from http://astronomy.swin.edu.au/~pbourke/other/interpolation/ */
static inline float CosineInterpolate(float y1, float y2, float mu)
{
   float mu2;

   mu2 = (1-cosine(mu*M_PI))/2;
   return(y1*(1-mu2)+y2*mu2);
}

/* source : http://astronomy.swin.edu.au/~pbourke/other/interpolation/ */
static inline float CubicInterpolate(float y0, float y1, float y2, float y3, float mu)
{
   float a0,a1,a2,a3,mu2;

   mu2 = mu*mu;
   a0 = y3 - y2 - y0 + y1;
   a1 = y0 - y1 - a0;
   a2 = y2 - y0;
   a3 = y1;

   return(a0*mu*mu2+a1*mu2+a2*mu+a3);
}


#endif

